Bistable Chimera Attractors on a Triangular Network of Oscillator Populations 
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We study a triangular network of three populations of coupled phase oscillators with identical 
frequencies. The populations interact nonlocally, in the sense that all oscillators are coupled to 
one another, but more weakly to those in neighboring populations than to those in their own 
population. This triangular network is the simplest discretization of a continuous ring of oscillators. 
Yet it displays an unexpectedly different behavior: in contrast to the lone stable chimera observed 
in continuous rings of oscillators, we find that this system exhibits two coexisting stable chimeras. 
Both chimeras are, as usual, born through a saddle node bifurcation. As the coupling becomes 
increasingly local in nature they lose stability through a Hopf bifurcation, giving rise to breathing 
chimeras, which in turn get destroyed through a homoclinic bifurcation. Remarkably, one of the 
chimeras reemerges by a reversal of this scenario as we further increase the locality of the coupling, 
until it is annihilated through another saddle node bifurcation. 
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I. INTRODUCTION 



While studying a continuum of identical oscillators on 



a ring with nonlocal coupling, Kuramoto et al. 13 1 
discovered a remarkable state where the population of 
oscillators splits into two subpopulations, where one is 
synchronized and the other is desynchronized. This 
state was later dubbed a chimera, alluding to a monster 
in Greek mythology that consists of incongruous parts. 
Since then, several groups have explored the nonlinear 
dynamics of chimera states [Ml M, [u|, [H QJ, H, H . 
Their emergence on the ring was first investigated ana- 
lytically by Abrams and Strogatz 0, S| ■ They found that 
chimera states were born through a saddle node bifur- 
cation, which appears to be the typical scenario for the 
emergence of chimeras on all network topologies inves- 
tigated so far. Shima and Kuramoto J33| showed that 
chimeras also exist on 2D lattices with free boundaries 
in the shape of spiral waves: here, the center of the spi- 
ral, characterized by a topological defect, is replaced by 
a desynchronized core with finite positive radius 23] . 

A recent breakthrough in the field of coupled oscilla- 
tor systems was made by Ott and Antonsen who showed 
that the dynamics of infinitely many oscillators could be 
reduced to a system of finite ordinary differential equa- 
tions. The method has since then been generalized and 
put into a formal mathematical context by Pikovsky and 
Rosenblum [3l[ and Marvel et al. j2(| independently, 
who show that the dynamics of each population may be 
reduced to a flow described by three variables plus con- 
stants of motion. Abrams et al. [l[ employed the method 
to investigate a system of two interacting populations of 
phase oscillators; they were able to calculate the saddle 



'Electronic address: erik.martens@ds.mpg.de 



node bifurcation for this system analytically and showed 
that the chimera states undergo a further change of sta- 
bility to become breathing chimeras via a supercritical 
Hopf bifurcation. 

Whereas many of these studies have been focusing on 
the nature of the chimera in idealized settings, others 
have been studying how chimeras would occur or behave 
in systems closer to real world situations: Omel'chenko 
et al. 28] show that a network of globally coupled os- 
cillators, subjected to delayed feedback stimulation with 
spatially decaying profile, also induces chimera states; 
thereby they argue that chimera states indeed are a 
generic feature of coupled oscillator systems. Delay cou- 
pled systems are also investigated by Sethia et al. [32| 
who discover the clustering of chimera states. Makovet- 
skiy et al. [2(| mention they found chimera-like states 
in systems of three level cellular automata related to las- 
ing systems in combination with spiral waves. Laing in- 
vestigates the robustness of chimera states against the in- 
troduction of nonidentical oscillators with heterogeneous 
frequencies in several systems For the case of two 

nonlocally coupled populations studied by Abrams et 
al. Laing shows how chimeras may both be desta- 
bilized and stabilized as the strength of heterogeneity 
(i.e. the width of the frequency distribution) of the os- 
cillators is varied [TH- Finally, we mention that in neu- 
roscience, spatially localized "bumps" of neural activity 
are found in networks of spiking neurons - such states 
have been proposed as mechanisms for visual orientation 
tuning and working memory, and have been related to 
chimeras [13, El- 

Important questions still remain. In particular, which 
topologies allow for chimera states? In other terms, 
can we classify the network structures that allow for 
chimeras? Even on one-dimensional domains the situ- 
ation is unclear. It has been shown that chimeras may 
exist on a ring [l3T | with a continuum of oscillators, but do 



chimeras also exist on a finite line segment or the infinite 
line, i.e. on a non-periodic domain? The observation 
that chimeras exist in the shape of spiral waves on 2D 
domains [IH [33| is a hint that this might be possible. 

One may seek to answer this particular question by 
lowering the dimensionality of the problem and discrctiz- 
ing a one-dimensional domain using oscillators popula- 
tions. As in Abrams et al. [lj|, one assumes that oscilla- 
tors within a given population are coupled more strongly 
to each other than they are to those in neighboring pop- 
ulations, thereby defining a spatial structure on the net- 
work. Thus the simplest network that exhibits both 
chain-like and ring-like character consists of three popu- 
lations, as shown in Fig. [TJ By tuning a new structural 
parameter c, a rotationally symmetric (ring-like or "tri- 
angular") network (Fig. Q] (b)) is deformed into a less 
symmetric, chain-like structure (Fig. [T] (a)). The im- 
pact of changing the network topology in such ways is 
discussed in a companion article [21] . 

This paper is wholly devoted to the study of the purely 
triangular network as shown in Fig. [T] (b). While exam- 
ining this simple case, we surprisingly found the coexis- 
tence of two stable chimera attractors. These attractors 
differ in the number of desynchronized populations (one 
or two). Their occurrence is unexpected because the ring 
with a continuum of oscillators 0,0], sharing the same ro- 
tational symmetry, has only a single stable chimera state. 
The finding of multiple coexisting chimera attractors is 
novel and is relevant in various contexts. First, if we let 
the number of oscillator populations go to infinity, we 
would retrieve the system of a ring studied by Abrams et 
al.; however, this system is not known to yield any multi- 
stable chimeras, which forms an interesting discrepancy. 
Second, the finding of bistable chimeras is noteworthy 
in the context of the question raised just recently [27| 
about what the basins of attraction leading to chimera 
and globally synchronized states actually are. Third, in 
connection with the above mentioned 'bump states', mul- 
tistability of chimera states may be relevant in the study 
of mathematical neural models, where so called switching 
processes 14 and the competition of attractors in working 
memory [9[ are of interest (see Discussion) . 

The organization of the article is as follows. We intro- 
duce the governing equations in Section II and summarize 
the derivation to obtain the reduced set of equations, as 
described in [l| . Next we derive the equations implied by 
special symmetries consistent with chimera states. These 
are analyzed in Section III, together with all their possi- 
ble bifurcations. Section IV discusses our results in the 
context of numerical simulations. Finally, our findings 
are discussed in Section V. Additional results on stability 
of some simple symmetric states for arbitrary networks 
are outlined in the Appendix. 



II. GOVERNING EQUATIONS 

The governing equations are given by 

3 N cr> 

<t' = 1 a j=l 

where the phases of the oscillators are defined by 6, i 
denotes the individual oscillators belonging to the pop- 
ulation with indices a = 1,2,3, each of which has N a 
oscillators, and where a is the phase lag parameter. 
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FIG. 1: Networks of three populations of oscillators. The gray 
disks symbolize the populations of oscillators, populated by 
individual oscillators symbolized by black dots. Their bidirec- 
tional coupling is represented by black arrows, (a) Chain-like 
general case with structural tuning parameter c. (b) Trian- 
gular network structure, corresponding to c = 1. Each pop- 
ulation has a self-coupling of unit strength 1, and is coupled 
to the neighboring populations with strength 1 — A. 

The coupling kernel K aa i describes the strength be- 
tween populations a and a' . The coupling strength is as- 
sumed to decay with increasing separation between the 
populations on the network. Within a population, the 
oscillators interact with strength K aa i — 1. Neighboring 
populations couple more weakly, with strength 1 — A as 
displayed in Fig. Q] (b) . We then have 

(l 1-A 1-A\ 
K aa , = 1 - A 1 1-4 . (2) 
\1-A 1-A 1 J 

In the case of A = 0, we retrieve the case of a glob- 
ally coupled network. Thus A quantifies how 'far' we are 
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from global coupling. This network has the same rota- 
tional symmetry as a continuum of oscillators on a ring, 
studied by @, y|, [l3| , and generalizes the problem with 
two populations discussed by Abrams et al. 



A. Reduction to low-dimensional system 

We consider the limit of infinitely large populations, 
i.e. N a — > oo. This allows us to reduce the problem to 
a finite set of equations, using the ansatz introduced by 
Ott and Antonsen [2{|, outlined below. In this limit, it is 
natural to describe the dynamics of the system in terms 
of the oscillator density distribution f a (9), which evolves 
according to the continuity equation 



df d 

dt de u ' 

The velocity of the oscillators is then given by 

3 r 27T 



(3) 



a)f (9',t)d9'. 







(4) 



To keep the notation simple we denote 9 a by 6 and 9 a i 
by 9'. It proves convenient to define a complex order 
parameter 
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e i6 'f(9',t)d9', 



(5) 



which defines a weighted average over all oscillators; we 
therefore refer to this as the global order parameter. Fol- 
lowing Kuramoto's footsteps [l2T - fl4l ]. we rewrite the ve- 
locity in terms of the order parameter and find 



v a = U) + Im 
I 



'f(9',t)d9' 



-id _ — ia 



2i 



i a 

e t0 e ta zt 
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(6) 



Following Ott and Antonsen [29( , we now restrict atten- 
tion to a special class of density functions in the form of 
a Poisson kernel 



5>„(t)e")' 



c.c. 



(7) 



where c.c. is the complex conjugate of the expression un- 
der the sum. The implications of this special ansatz and 
its validity will be explained in more detail in Section fVl 
Substitution of / CT and v a into the continuity equation ([3]) 
yields an exact solution, so long as a a evolves according 
to 







a a + iLoa a + -a a z a e 



(8) 



It remains to express the order parameter in terms of this 
ansatz. We find 



r (f) - ^X CTCT <,(i). 



(9) 



<t' = 1 



Finally, we express the amplitude a a in polar coordinates 
as 



a a = p a e 



-%<t>a 



By division of (HJ by e , we obtain 



= p a - i$oPo + iupa 

1 3 

+ 2 Pa Zl K °°'P°' 
cr'=l 

1 3 

~o H K aa ,p a ,e 



i{<fr a >~<t><T-a) 



and by separation of real and imaginary parts we find 

1 2 " 

1 -Pa 



Pa 



^ Kea'Pa, Sm {<j> a i - CT + 0), 



cr'=l 
I + p 2 3 

4>cj = W Kaa'Pa' COS ((/v ~ 0<r 



where we introduce the definition 



/? = tt/2- a. 



(10) 



(11) 



These equations describe the dynamics of our system in 
terms of the variables a a . Notice that in contrast to 
z a , they do not represent averages over all populations, 
and therefore, we refer to them as local order parame- 
ters. Thus any synchronized population a of oscillators 
is characterized by p a = 1. (These results are trivially 
generalized to the case of a network with arbitrarily many 
populations a = 1,2, ... , N.) 

In what follows, we will make particular assumptions 
about the symmetries of the solutions that we expect 
to find. Then we will analyze existence, stability and 
bifurcations of the states of interest. 



B. Manifold of symmetric states (SDS and DSD) 

To make progress on analyzing the chimera states, we 
will have to make certain symmetry assumptions. Per- 
fectly synchronized populations have p a — 1; desynchro- 
nized populations have p a < 1, and consist of oscillators 
that drift relative to one another and to the synchro- 
nized populations. Let S and D denote synchronized 
and desynchronized populations, respectively. Then in a 
triangular network, we can distinguish only two chimera 
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states, namely SDS (sync-drift-sync) and DSD (drift- 
sync-drift); all other permutations of S and D give equiv- 
alent states because of the rotational invariance inherent 
to the triangular network. 

Another class of solutions can be described as SSS, 
corresponding to a globally synchronized state where all 
oscillators are in sync, but with different synchronized 
mean phases 0;. We may distinguish three cases, i.e. 
0i = 02 = 03, 4>i — 4>3 02 and the state where all 
phase angles are different. These solutions are analyzed 
in Appendix A. These states are of less interest to us and 
we restrict ourselves to chimera states and their emer- 
gence in parameter space in the following. 

The symmetry of our coupling kernel ([2]) in combina- 
tion with pi = p3 implies that 0i = 03 and hence popu- 
lations 1 and 3 are phase-locked. The SDS state is then 
defined via 

Pi = P3 = 1 and p = p 2 < 1, 

01 = 03, 

whereas the DSD state is given by 

p = pi = p 3 < 1 and p 2 = 1. 

01 = 03, 

We define the phase difference of the angular order pa- 
rameter between the synchronized and desynchronizcd 
states by 

-0 = 01 - 02 = 03 - 02- (12) 

Applying these symmetry assumptions to (|10p and sub- 
stituting the coupling kernel defined in ^ , we obtain the 
equations describing the SDS states 

1 — 2 

p = -^_-[2(l-A)sin^ + /3)+psin/3], 

ip = -(2-A)cos / 0-(l-A)pcos(-^ + y S) 
1 + p 2 

+ 2(1 -A) cos (ip + /3 +pcos/3 , (13) 

2p 

and the DSD states 

p = ±^[(2-A)psmP + (l-A)sm(-1> + P)], 
1 + P 2 

i) = — [(2- A)p cos /3 + (I -A) cost -ip + /3)] 

2p 

+ 2(1 - A)pcos{4) + P) +cos/3. (14) 

Note that these equations hold only if we restrict atten- 
tion to symmetry-preserving perturbations. The fixed 
points of (|13ll4p correspond to phase-locked solutions of 
the original system (at the macroscopic level of the local 
order parameters). 

By reduction of the full system of oscillators (fTJ) to a 
low dimensional system for the local order parameters, 
we have cast our problem into a two dimensional system 
represented by Eqs. (|13ll4j) . This enables us to study 
the problem in the phase plane. 



III. ANALYSIS 
A. Phase portraits 

Unfortunately, the equations for the SDS and DSD 
states cannot be solved in closed form. Before we get 
deeper into the matter of analyzing these states, let us 
get a quick intuition of their behavior by inspecting the 
phase planes of their corresponding equations, which will 
guide us in the subsequent analysis. Their phase por- 
traits shown in Fig. [5] and Fig. 0] represent a sweep in 
parameter space with increasing values of A while keep- 
ing the value of /? = 0.1 constant . 

Consider first the SDS symmetry. For small values of 
A (close to global coupling) , we only observe a nodal sink 
and source on the unit circle; these points correspond to 
the in-phase SSS solutions (Fig. Ufa)). Increasing A 
further, a saddle-node pair is born very close to the unit 
circle. For larger A, the node moves closer to the origin, 
implying that the order of the desynchronized popula- 
tion decreases and the chimera state becomes more pro- 
nounced (if we instead increase the values of /3, this crit- 
ical point starts to move clockwise while getting closer 
to the origin). The node has become a stable spiral 
(Fig. Hlb)), and at a critical value of A, it loses stabil- 
ity through a Hopf bifurcation and a limit cycle is born 
(Fig. [He)). The amplitude of the order parameter p of 
the drifting population starts to oscillate, and is there- 
fore called a breathing chimera. As we raise the value of 
A more, the limit cycle gains in amplitude until it col- 
lides with the saddle: the limit cycle is destroyed in a 
homoclinic bifurcation (Fig. [H[d)). The resulting bifur- 
cation diagram is shown in Fig. [3] (a). The saddle-node, 
Hopf and homoclinic bifurcation curves all intersect in a 
Takens-Bogdanov point with codimension two. 

For the DSD symmetry, we observe a scenario that is 
qualitatively similar to the previous case. However, sur- 
prisingly, the whole scenario appears a second time in the 
upper part of the parameter plane, but now in reversed 
order, as shown in Fig. [3](b). We again increase A while 
keeping (3 constant, as shown in Fig. [4] For all param- 
eter values, we find two synchronized SSS solutions on 
the unit circle: one is a nodal sink in (p,ip) = (1,0), 
but what in the SDS case before was a nodal source, 
is now a saddle. Also notice that a new fixed point has 
appeared in the left half of the unit circle in the form 
of a spiral source (Fig. Ufa)). This is the second, cur- 
rently unstable, DSD chimera seen in the upper half of 
the bifurcation diagram. As A increases and the cou- 
pling becomes more local, a saddle-node pair is born in 
the right half of the unit circle (Fig. 0|b)); again, its node 
then becomes a spiral and loses stability as the coupling 
strength becomes more local, and the resulting limit cycle 
(Fig. HJc-d)) gets ultimately destroyed in a homoclinic 
bifurcation (Fig. 0Je)). Whereas one chimera has been 
rendered unstable, we observe that, above a critical A, a 
stable limit cycle has formed around the spiral source on 
the left half of the circle: the twin of the DSD chimera in 




FIG. 2: Phase portraits for the SDS chimera, with increasing values of A at constant /3. Real (a;) and imaginary (y) components 
of a( CT ) are shown. The unit circle is displayed in gray. Stable fixed points are shown as solid (red) and half-filled (green) circles, 
respectively. Limit cycles are emphasized in blue color. Stable and unstable manifolds are shown as red solid and green dashed 
trajectories, respectively. The point in (p, ?/>) = (1, 0) is a nodal sink. The position of the nodal source depends on ft and moves 
in clockwise direction with growing values of j3. 



its breathing mode has emerged (Fig. IUc)). From here, 
the bifurcations happen in reversed order, the source be- 
comes first a spiral node (Fig. 0Jf)), i.e. a stable chimera, 
which is eventually annihilated in a saddle-node bifurca- 
tion. The resulting bifurcation diagram is seen in Fig. [3] 
(b). 

B. Calculation of bifurcation curves 

In order to calculate the saddle-node and Hopf curves 
of the SDS solutions, we must linearize (|13[) around the 
appropriate fixed point. This task amounts to solving 
the fixed point equations implied by Eqs. ([TU]) and (fTl)) 
simultaneously with the saddle node condition, 

det(J) = 0, 



or with the Hopf condition, 

tr(J) = and det (J) > 0, 

where J denotes the Jacobian of (| 13|) or , respectively. 
For the SDS symmetry, we have 

Jn = \ [(l-3 /9 2 )sin/3-4(l-^)psin(/3 + V)] , 
J12 = (l-A)(l-p 2 )cos(/3 + ^), 
J21 = -o [cos /3(p 3 - (1 -A) cos ip) 

- sin/3sini/;(2p 2 - 1)(1 — A)] , 
A — 1 

J22 = [(1 + 2p 2 ) cos i/)sin/3 + cos /3 sin • 
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(a) 



A 



breathing chimera 
stable chimera 




FIG. 3: Bifurcation diagram for the SDS (a) and DSD (b) 
chimera. The displayed curves are: the saddle-node curve 
(solid red), the Hopf curve (dashed blue), and the homoclinic 
curve (dotted black). Dots mark the bifurcation points ob- 
tained by inspection of the phase plane. The homoclinc curve 
(dotted) is an interpolation based on these points, whereas all 
the solid curves were obtained analytically. 



The fixed point condition implied by (fT3j) yields the non- 
trivial solution 



p = 2(A- l)csc/3sin(/3 + ip). 



(15) 



Substitution of this expression into the fixed point con- 
dition for i/j results in 

= (A - 2) cos /3 + i esc (J3 + i/j) sin ip 

+ (1 - A) 2 esc P [sin 2/3 + 2 (cot /3 sin 2 ip + sin 2ip^L6) 

Unfortunately, this equation cannot be solved in closed 
form for if>, which in turn would allow us to express A in 
terms of /3. We settle therefore for a series approach in 
A and ip, as follows: 



JV 



.4 



N 

]T a k (3 k 



■0((3 N+1 ) 
-0((3 N+1 ) 



(17) 
(18) 



We substitute these two expressions into fixed point equa- 
tion (TTB1) and the saddle node condition, and solve the 
resulting equations for each power of ft. This leads to 
the following expression for the saddle node curve: 



35283 

"80 
2949379 



/3 5 



210247 



160 



■0 8 



2355 

16 P 
872617 



224 



256 



2744116261 
80640 



0(/3 10 ). 
(19) 



Using the Hopf condition and proceeding in the same way 
we also find the Hopf curve, approximated by 



A H (J3) = 0.447153+ 1.34639 /3 2 



3.34 371 /3 4 + C>(/3 6 ). 

(20) 



The Takens-Bogdanov point is determined by numeri- 
cally solving the saddle-node, Hopf and fixed point con- 
ditions simultaneously. It is located at 



(0,A)sds = (0.1974,0.5092). 



(21) 



The bifurcation curves for the DSD chimeras are ob- 
tained in an analogous procedure. Solving the equations 
expanded in series is now a bit trickier due to the co- 
existence of two branches. For brevity, we shall only 
summarize our findings. The two saddle node curves are 
approximated by 



.4 
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08) 



9/3-72 /3 2 
1130943 



1059 



/3 3 - 38 55 (3 A 
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/3 5 

40 

854234093 



1039276 



P 7 



P 6 
o 

78311783 



560 

3309788681161 



.4 



SN,2 



08) 



40320 



7 

- OG3 10 
6421 



/3 8 



120 



(22) 

0(/3 6 ). 
(23) 



Unfortunately, the series approach for the lower second 
branch converges extremely slowly and doesn't match all 
the way to the Takens-Bogdanov point (even going to this 
high order doesn't help!). However the effort to find the 
series coefficients was not all in vain: a Pade approximant 
based on the above power series at order three does an 
excellent job in matching the data points retrieved from 
the examination of the phase portraits. We have 



A 



SN,2 



08) 



9 /3 



i , 2580226 a 
^ 301403 V 



1521018 a2 
301403 P 



90123833 02 
6028060 P 



, 287446827 o3 
' 6028060 P 



200207828 03 ' 
4521045 P 

(24) 



fc=0 



Finally, the Hopf curves are approximated by 
A H ,i(P) = 0.593737+ 1.14491 (3 2 



7 




FIG. 4: Phase portraits for the DSD chimeras, with increasing values of A at constant ft. Real (a;) and imaginary {y) 
components of a( CT ) are shown. The unit circle is displayed in gray (the unstable green dashed manifold of the saddle coincides 
with it). Stable and unstable fixed points are shown as solid (red) and half-filled (green) circles, respectively. Limit cycles are 
emphasized in blue color. Stable and unstable manifolds are shown as red solid and green dashed trajectories, respectively. 
The point in (p, ip) = (1,0) is a nodal sink. The position of the saddle depends on /? and moves in clockwise direction with 
growing values of ft. 



+4.55308/3 4 + O(/3 6 ), (25) 
AhAI 3 ) = 0.885408- 1.15074 /3 2 

-3.96289 f3 4 + 0{f3 e ). (26) 

The Takens-Bogdanov points are located at 

(P,A) d bd,i = (0.2132,0.6615) (27) 
(/?, A) DSD , 3 = (0.1903,0.8359). (28) 

For all symmetries, we found an excellent agreement of 
our perturbative results with the bifurcation points ob- 
tained from inspection of the phase portraits. We did 
not derive an analytical expression for the homoclinic bi- 
furcation curves; the curves shown in the Figs. [3] (a) and 
[3] (b) are based on data points obtained from inspecting 
the phase portraits while varying the parameters. 



IV. NUMERICAL SIMULATIONS 



We have obtained analytical results describing the dy- 
namics of our triangular network of oscillator popula- 
tions, using two different reductions: firstly, we reduced 
the governing Eqs. (Q]) using the Ott Antonsen method. 
And secondly, we have assumed that the system attains 
certain symmetry states that allow for chimera states; 
these symmetries need not be transversely stable to per- 
turbations off the symmetry manifold. Thus, the equa- 
tions obtained from these reductions may not necessar- 
ily account for the complete dynamics of the governing 
equations, and we checked if our analytical results agree 
with numerical simulations of the governing equations. 
We did this with a finite, but what may be considered a 
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sufficiently large oscillator population (N^~ 40 ). 

We used two different methods to generate initial con- 
ditions that would lead to the appearance of chimera 
states: Firstly, for the phases of the desynchronized pop- 
ulations we used an initial condition in the shape of a 
bump, specifically, a Gaussian distribution in the shape 
of 4>1' d ~ exp (—j(i/No- — 1/2) 2 ) with an appropriately 
chosen decay rate 7. The second method was chosen 
with the intention to place the system right on the Ott- 
Antonsen (OA) manifold. This was accomplished by gen- 
erating a phase distribution that is consistent with the 
Poisson kernel ©. For both methods, the phases for 
the synchronized populations are given by (j\ ,s — 0. To 
achieve this, we solved the SDS and DSD Eqs. (TT3|) and 
(1141) . respectively, for fixed points (p,ip). This in turn, 
enables us to compute the Poisson distribution f a (9,t) 
defined by (J7J), by using the definition of the order pa- 
rameter, a a — p a e~ 1 ^" . Because this function defines 
the probability with which oscillators populate a certain 
phase, we may use its inverse cumulative distribution 
function to construct from it a set of phases that is con- 
sistent with the OA-manifold. The system should remain 
close to the OA-manifold, because of its invariance. Un- 
less mentioned otherwise, we used this latter method. 

We first confirmed that the unreduced system would 
exhibit all types of chimeras predicted by our analysis, 
and that they would correspond to stable states, for var- 
ious points in parameter space. These states were ob- 
served with both N a = 20 and N a = 200 oscillators per 
population. The observed behavior is that the system 
first goes through a tiny transient and reach an attract- 
ing state which was confirmed to be stable even for long 
computation times. The transient may be explained by 
the fact that the system, due to its finite size, can only 
be approximately on to the OA-manifold (or a member 
of the OA- family, as explained in the Discussion). 

Next we checked if the critical parameter values of sad- 
dle node, Hopf and homoclinic bifurcation in the full sys- 
tem would be in accord with the critical values obtained 
from our theory, see Figs. [5] (a) and (b). In order to do 
so, we held the value of j3 fixed and continued a solution 
through 20 increasing values of A. To initiate the con- 
tinuation we used an initial condition consistent with the 
OA-manifold. 

In Fig. [5J we show fixed points and oscillation am- 
plitudes of the order parameter, obtained analytically 
from (IT3l) and (IT4l) . as dashed curves. Our simulation 
results are superposed, and were obtained as follows: To 
remove transient effects from our analysis, we only con- 
sidered the last 2/5 of the computed time series. Instead 
of only detecting the global maximum and minimum of 
the series, we detected local maxima, shown as light gray 
dots, and minima, shown as dark gray dots. We chose 
to do this, because it would allows us to see the traces 
of new appearing periods, that potentially may occur in 
the unreduced, highly dimensional system. However, fi- 
nite size fluctuations also cause small oscillations; this is 
the reason we see some small amount of blurring in the 



data (For similar reasons, maxima and minima on the 
stable branch do not coincide.) 
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FIG. 5: Bifurcation diagram obtained from numerical simu- 
lation, for states SDS, shown in (a), and DSD in (b). Light 
dots and dark circles correspond to local maxima and minima, 
respectively, that are detected by the algorithm. The dashed 
curve represents the analytical result for the continuum case 
(N 00). The computations were performed with N a = 40 
oscillators per population for a simulation time of T — 100. 
The question mark indicates that we could not conclusively 
confirm the existence of the second DSD state. (The kink 
in the lower left dashed branch is an artefact from the limit 
cycle reaching into the left hand side quadrants, as p is not 
measured relative to the limit cycle center but to the origin.) 

Despite these undesired effects, we can clearly demon- 
strate the onset of the Hopf and the homoclinic bifurca- 
tions in the simulation. Consider first the SDS chimera, 
shown in Fig. [5] (a). We expect that finite size effects af- 
fect the locations of all the bifurcation points. While we 
would have to compute at a higher resolution for the con- 
tinuation to actually see this happen for the Hopf bifur- 
cation, it is more apparent for the homoclinic bifurcation: 
the limit cycle oscillation brings the system periodically 
close to the saddle point (as seen in the phase plane); 
therefore, the larger finite size fluctuations are, the more 
likely the system is to be kicked off the limit cycle and 
is instead attracted to the nearby all-in-phase SSS state 
on the unit circle. This is seen in Fig. [5] (a) , where the 
limit cycle oscillation disappears at a smaller value of A 
than it does for the analytic result (sold curve indicates 
SSS state). The same behavior is observed for the DSD 
chimera in Fig. [5] (b), but much more pronounced: the 
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system snaps to the in-phase state (solid curve) much ear- 
lier than expected for a continuous system of oscillators. 
Similarly, we can continue the solutions in reverse direc- 
tion, and check that the chimera states are annihilated 
via a saddle node bifurcation. 

Whereas we easily managed to show DSD states in 
the lower half of the parameter plane, our attempts to 
place the system on the second DSD attractor had little 
success: for all trials the system would eventually reach 
the in-phase state. Increasing the number of populations 
didn't seem to remedy the matter, as one might be led 
to think from our previous experience; rather, the sys- 
tem would stay obediently in the DSD state for a little 
while and then suddenly jump into the in-phase state. 
However, it doesn't seem to be entirely clear whether the 
second attractor of the DSD state is inherently unsta- 
ble, or if it is just very hard to stay on that manifold due 
to a combination of factors, given by finite size effects 
and the location and shape of the attractor (see phase 
plane in Fig. Specifically, in the case of a breathing 
state, the limit cycle is always very close to the invariant 
field defined by the in-phase state on the unit circle. In 
conclusion, it is likely that this second DSD attractor is 
unstable to symmetry breaking perturbations. 



DISCUSSION 



We have investigated a triangular network of popula- 
tions of sinusoidally, nonlocally coupled oscillators with 
identical frequencies. Our analysis approximates the 
practical case of large populations by considering the 
limit of infinitely many oscillators per population, for 
which we may use the Ott-Antonsen method to reduce 
the governing to a finite set of equations [2{|. By as- 
suming symmetries compliant with chimera states, we 
further reduced these equations and studied their emer- 
gence in a phase plane analysis. Saddle node and Hopf bi- 
furcation curves are determined using perturbation tech- 
niques, and homoclinic bifurcation curves by observation 
of phase portraits. The resulting stability diagram is a 
variation of the one found by Abrams et al. but 
is governed by the emergence of two different types of 
chimeras: one with a desynchronized 'zone' of narrow 
width (SDS), and the other with a desynchronized zone 
of twice the width (DSD). DSD chimeras exist in two 
regions of parameter space, and differ in their mean (dif- 
ference) phase angle ip; for the DSD state near global 
coupling (small A), ip remains close to zero, and for the 
other near local coupling (large A) closer to tt. Numerical 
experiments demonstrate that chimeras persist for small 
oscillator populations, and that the SDS and the DSD 
chimera near global coupling are truly stable attractors in 
the unreduced system described by Eqs. ([TJ. Finally, we 
observe for the first time the feature of bistable chimeras. 



A. Reduced system and numerical experiments 

Ott and Antonsen [30] showed that the long time dy- 
namics of a large class of phase oscillator systems (includ- 
ing this one) is attracted to the manifold defined by (|7|) if 
the natural frequencies are drawn from peaked distribu- 
tions. This result has been confirmed in several studies 
H 0, EH HIJ . The situation is somewhat different if we 
consider the case of identical frequencies. Pikovsky and 
Rosenblum [3l[ , and Marvel et al. 24 1 have shown inde- 
pendently that the dynamics of each population can be 
reduced to a flow described by three variables plus con- 
stants of motion. The unifying picture is that for identi- 
cal frequencies one has a whole one-parameter family of 
invariant manifolds, including the OA manifold. These 
manifolds are neutrally stable with regards to perturba- 
tions in transverse direction to themselves (such pertur- 
bation would result in placing the dynamics into a neigh- 
boring manifold). Once the frequencies are spread this 
family collapses and the OA manifold is left. 

Because this degeneracy in the case of identical fre- 
quencies might lead to unanticipated results, and because 
the symmetry manifolds defined by SDS and DSD don't 
need to be stable to perturbations in transverse direc- 
tions, we checked that the chimera states are true attrac- 
tors in the unreduced system, Eqs. (JTj), by numerical sim- 
ulation. The sequence of bifurcations with saddle node, 
Hopf and homoclinic bifurcations is indeed reproduced in 
the unreduced system and the bifurcations appear close 
to the predicted critical values, even though the number 
of oscillators per population is small. The homoclinic bi- 
furcation makes an exception in the sense that it seems 
to happen already for small A; we have argued that this 
likely is an artefact for large limit cycles, where the tra- 
jectory gets kicked off the orbit because of finite size fluc- 
tuations in combination with the appearance of a nearby 
saddle point. All chimera states, with exception of the 
DSD state near local coupling, have proved to be true 
attractors in the unreduced system. 



Relation to heterogeneous frequency 
distributions 



A recent study by Laing [15| generalizes the problem 
of a network with two oscillator populations investigated 
by [H to the case of heterogeneous frequencies. Laing 
showed that for this and various other network topolo- 
gies, the chimera state is robust - within limits - to het- 
erogeneity in the intrinsic frequencies of the oscillators. 
In particular, he finds that the chimera state remains sta- 
ble for populations with nearly identical oscillators, that 
is, with a narrow width of the distribution. The stability 
results obtained from our analysis should therefore be the 
same as the one obtained for the dynamics of oscillators 
with almost identical frequencies. 
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C. Aperiodicity and chaos 

In our setting of the problem with identical oscillators, 
we were able to find solutions that do not lie on the OA- 
manifold (to find them, use initial conditions lying off 
the OA-manifold, e.g. Gaussian 'bumps' or add noise to 
the OA-compliant initial conditions). Such states may 
be related to the quasi-periodic states shown to exist by 
Pikovsky et al. [3l|, but perhaps also to chaotic states. 
We haven't gone further into this question, but have no- 
ticed irregular behavior in some of our simulations that 
look aperiodic, and we think it might be possible to find 
something like a chaotic chimera state. 



D. All-in-phase states 

In the case where all populations are fully synchro- 
nized, i.e. p a = 1, the dynamics of the system effectively 
becomes the one governed by three coupled oscillators. 
The stability of such states is analyzed in the Appendix. 
This case has already been studied in the context of three 
coupled sinusoidal limit cycle oscillators by Mendelowitz 
et al. [25j and for relaxation limit cycle oscillators by 
Bridge et al. Q, who found that two rotating waves, 
clockwise and counter-clockwise rotation, are possible. 



E. Relation to ring of oscillators and network with 
two oscillator populations 

We mention the connection of our study to the work 
done by Abrams et al. [l| , who studied a similar system, 
but with only two oscillator populations. With its two in- 
phase-locked populations, it is unclear whether the SDS 
or the DSD state compares best to their - using our 
terminology - SD chimera. Certainly, it can be said they 
both act like a two population system in disguise. In this 
sense, the two population system is a degenerate case of 
our triangular network. 

The same authors also studied a continuum of oscil- 
lators on a ring 0, Q. In our picture with oscillator 
populations, this system is approximated by an infinite 
set of populations arranged on a ring. For the continu- 
ous ring, only a single chimera is known. Both SDS and 
DSD chimeras effectively act like a system made of two 
populations, but differ in their 'width' of drifting pop- 
ulations. In the case of the ring, this width is slightly 
larger than the synchronized region; from this point of 
view, the DSD chimera seems to match their state best. 

One may ask how the behavior of discrete ring-like sys- 
tems is affected as we increase the number of populations. 
One could imagine that more and more chimera states get 
added as we increase the number of populations, consti- 
tuting competing multistable attractors. But then, what 
happens as we take this continuum limit? Do they dis- 
appear, collapse such that only one of them, specifically, 



the one discussed in |3j, wins the competition, and dom- 
inates over all others and remains stable? In the case 
of the triangular network, the SDS and the first DSD 
chimera are truly stable attractors in the unreduced sys- 
tem, and so the situation stays inconclusive. The creation 
and annihilation of chimera states with varying numbers 
of oscillator populations is also of importance if we want 
to characterize the basins of attractions in complex oscil- 
lator networks, an issue pointed out recently by Motter 



F. Multistability and neural networks 

Our numerical experiments show that the existence of 
chimeras not only exist in the thermodynamic limit of 
N — > oo, but also for relatively low numbers of oscillators 
per population of N ~ 20. In view of the robustness of 



chimera states towards heterogeneity [15| , one may there- 
fore expect that chimeras also may emerge in biological 
settings such as neuronal networks. In this case, the co- 
existence of chimera attractors would have implications 
for the study of neuronal networks beyond their similarity 
to spatially localized " bumps" of neuronal activity men- 
tioned earlier. For instance, one may envision that mul- 
tiple coexisting chimera attractors encode states of asso- 
ciative memory 0,H| or play a role in decision making, 
where different initial conditions lead to different states. 
The characterization of basins of attractions for chimeras 
therefore represents an interesting issue, and represents a 
topic in the field that is not well explored 27]. Alterna- 
tively to stable-state dynamics, some neuronal models, 
based on the mechanism of switching among unstable 
saddles, have drawn attention lately [4J. In these models, 
heteroclinic connections generate characteristic patterns 
of neural activity and thus represent information. One 
may conceive oscillator networks with the possibility of 
switching between chimera states while varying system 
parameters. Provided that such chimera sequences, en- 
coding different spiking activity patterns, are realizable, 
they may constitute an interesting field of research. We 
note that in this context oscillator populations, or rather 
neurons, represent redundancy in computing networks, 
which significantly enhances their reliability, as shown 
already by von Neumann [34j . In a recent experimental 
study Feinerman et al. Q build neuronal computational 
devices and demonstrate how reliability is significantly 
improved by employing neuronal ensembles for each logic 
unit, similar to the oscillator populations. 
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Appendix A: Stability of Symmetric States 
1. Existence of Symmetric States SSS 



The first eigenvalue is a manifestation of the rotational 
invariance of the system (the system only depends on 
phase differences and is effectively five dimensional). 



The states we consider here are fully synchronized, i.e. 
p a = 1 for all a = 1,2,3. We consider the symmetries 
where 0i — 03 . In order to derive a fixed point condition 
for these states, it is therefore sufficient to consider the 
equations that are specialized to the symmetries SDS 
and DSD given by Equations (jT3"]) and (|14")) . Applying 
the fixed point equation to either of them yields the con- 
dition 



We first compute the stability of the trivially symmet- 
ric state defined by ip = 0. This degenerate state has the 
eigenvalues 



= A(A — (2 A — 3) sin^) 3 (A - (A - I) sin/3) 2 (A7) 



= (1 — A)(cos(3 — cos/3cosi/; + 3 sintp sin/3), (Al) 

which is reduced to these cases: 

A = 1, (A2) 

P = with cos-0 = l, sin ip ^ 0, (A3) 

= 0, (A4) 

a 3sin/3sin?/; 

cosp = with cos-0 T 1- (A5) 

cos ip — 1 

The first two conditions are the degenerate cases repre- 
senting the A and j3 axis. The third condition corre- 
sponds to the fixed point (p,V0 = (1,0), and the last to 
the position of the node that is constrained to move on 
the unit circle in the phase portraits, see Figs. [5]and|31 

2. Stability of SSS States 

The computation of stability of these points is ac- 
complished by computation of the Jacobian of the six 
dimensional system (flOl) . with the coordinate system 
(pi, P2, P3, 0i, 02, 03), using a computer algebra system. 
Applying our symmetry assumptions, we find the eigen- 
values 

/ ° \ 

-sin^ + 2(A- l)sin (B + 0) 

(A - 2) sin/3 + (.4-1) sin (/?-</,) 
Ak ~ (A-2)sin/3 + (A-l)sin(/3-V0 ' [ ' 
(A-l)(2sin/3 + sin(/3-V)) 
\ (A - 1 ) ( 3 cos V> sin /3 + cos p shrtp ) j 



If we only consider parameter values A 6 (0, 1), a suffi- 
cient condition for linear stability is given by p £ (0, n). 

The less trivial state with vb ^ must be considered 
in combination with the fixed point condition (|A5I) . The 
signs of these eigenvalues were not obtained analytically 
but rather by graphing the maximal eigenvalue in the 
(/3, A)-plane. It turns out that this state has at least one 
positive eigenvalue except for the degenerate case where 
P = 0, 7r, and is thus always a saddle. This result is 
consistent with the behavior observed by inspection of 
the phase portraits of Eqs. (fT3")) and (|14p (in the case of 
SDS symmetry, the nontrivially symmetric state changes 
stability at P = w, but this holds only within the SDS- 
symmetry manifold, and has nothing to do with stability 
in the six (or five, for that matter) dimensional space.). 

It is possible to obtain a similar stability result for 
the general case of a network with arbitrarily many pop- 
ulations N, for the trivially symmetric point satisfying 
p CT = 1 and 0i = 02 = ... = 0jv, using Gershgorin's circle 
theorem. In this general setting, the point also becomes 
linearly unstable as /3 > n (provided that all row sums of 
the coupling matrix are strictly positive). Unfortunately, 
no similar result was obtained for the remaining N — 2 
fixed points that may occur on the unit-sphere related to 
the general problem. The calculation is tedious and will 
not be represented here. 
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